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Abstract 


In the present work mathematical formulation of /-integral and its finite element 
implementation is presented. A Matlab code is developed for evaluating /-integral for 
two dimensional crack problems and is extended to determine three dimensional point- 
wise /-integral along a curved shape crack front. The results from the code are verified 
against standard elastic problems for which the analytical results are known. 

Though the analysis of the components in nuclear industry involves geometric as 
well as material nonlinearity and the domain of /-integral is the elasto-plastic analysis, 
but as starting phase of the development of code, the /-integral is evaluated and validated 
in the context of Linear Elastic Fracture Mechanics (LEFM). 

A parametric study is carried out for a pipe having circumferential semi-elliptical 
surface crack for different aspect ratios in the elastic region. 
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CHAPTER 1 
INTRODUCTION 


1.1 GENERAL INTRODUCTION 

Per capita power consumption is a measure of living standard in today’s world. The 
major sources of power generation are fossil fuel, hydro power and nuclear power. While 
the reserves of fossil fuel are limited, the resources of the hydro power are highly 
dependent on the geographical topology, further, they are also limited. The share of 
nuclear energy is growing in power sector. In fact, it will be a major constituent of the 
energy sources in future. 

The safety is of utmost concern in design and manufacturing of components of 
nuclear power plant. But at the same time the cost of power generation must be 
affordable. Loss of coolant accident (LOCA) is most likely to occur in the pressurized 
water reactor (PWR) plant. The scenario referred to as “thermal shock” involves the 
possibility of fracture of a nuclear reactor pressure vessel during a loss of coolant 
accident (LOCA). 

One of the most likely causes of a loss of coolant accident in a nuclear power 
plant is a rupture in the piping system caused by stress corrosion. Many incidents of 
stress corrosion cracking have been reported. Consequently there is a great need for a 
quantitative understanding of the behavior of cracked pipe under normal operating and 
accident conditions — for example a seismic event. In most instances, the concern is for 
surface cracks that initiate at the inner surface of the pipe in the heat affected zone around 
a girth weld. Affected by the weld induced residual stress, these cracks tend to grow 
circumferentially and radially, sometimes attaining a size that is a significant fraction of 
pipe wall area. Fig. 1.1 shows an example of the cracks that were discovered in the wall 



of 4-in diameter stainless steel pipes in the boiling water reactor (BWR) in 1974. This 
example shows that substantial crack sizes can be achieved before detection. 




Figure l.lCrack detected in 4-in diameter recirculation by-pass of a boiling water reactor plant 

A prime consideration in analyzing cracked nuclear power plant pipes is to 
determine, if failure actually occurs, whether it will lead to a “leak-before-break- 
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condition”. In general the leak-before-break concept refers to a pressure containment 
system failure event in which a part-through-wall crack extends to become a through wall 
crack, thus allowing contained fluid to escape. If no further crack growth occurs, then the 
loss of fluid can presumably be detected in one way or another and the system can be 
shutdown safely. The alternative - the through-wall crack propagates along the wall - 
will lead to a catastrophic event. Obviously it must be avoided. If failure occurs, it is very 
desirable to assure that it will be confined to the leak-before-break mode. 

In nuclear power plant applications, it is necessary to show that the leak-before- 
break concept is the applicable failure mode in piping systems, where cracking has 
occurred or even could occur. Specially, it must be established that a pipe crack will be 
revealed by leak detection techniques before it reaches a condition where fracture could 
occur under normal operating or postulated accident condition. The design basis for 
accident used in nuclear plant regulation around the world is the so called full guillotine 
offset break (i.e., instantaneous circumferential fracture). This extreme condition has 
resulted in the incorporation of massive pipe whip restraints into nuclear piping systems. 
These restraints are not only very expensive to design and install, but they can reduce the 
reliability of in-service inspection while increasing the radiation hazard in the inspection 
process. The necessity of such devices in design stage and as modification in operating 
plants would be substantially relieved if leak-before-break conditions could be 
demonstrated. Consequently, there is currently a great deal of research interest in 
developing more precise fracture mechanics analysis methods. 

The material of nuclear reactor pressure vessel is low alloy steel while the 
material of piping is stainless steel. It is difficult to join the two different materials 
because of their different crystalline structure. For joining the feritic steel and austenitic 
steel a thin layer of a different material, which is weldable to both the materials, known 
as buttering material is applied to one of the materials and then welding is done. This 
joint is obviously having lesser strength than a joint between the same materials. The 
joint between two materials having different crystalline structure is known as bimetallic 
joint. So this bimetallic joint between the pressure vessel and the piping is more 
vulnerable than other joints in the remaining piping. 
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1.2 REVIEW OF LITERATURE 

The nuclear plant piping materials are very ductile and tough. In the materials selected 
for such service, crack growth is generally preceded by substantial crack-tip blunting 
while significant amount of stable crack growth can occur prior to the fracture instability. 
Linear elastic fracture mechanics techniques usually provide a very conservative 
prediction in such circumstances, conservatism that often prompts unnecessary remedial 
action. So J-integral is used to characterize the crack in the bimetallic welded pipe joint 

/-integral is the most used crack tip integral in fracture mechanics. Its role in 
nonlinear fracture mechanics was introduced by Rice [1, 2] in 1968, who provided the 
interpretation of / as a measure of the intensity of deformation at notch or crack tip and, 
in the context of inelasticity, as energy release rate. Within the thermodynamics 
framework, the integral was introduced by Cherepanov [3] in 1967 as an extension of 
Griffith’s approach to inelastic solids. Subsequently, the connection between / and 
Eshelby’s energy momentum tensor was noted [4, 5]. Translational and nontranslational 
conservation integrals have been presented by Knowles and Sternberg [6]; the energetic 
force interpretations were made by Budiansky and Rice [7]. In this context /-integral is 
member of translational conservation integrals. However /-integral possesses additional 
features which make it unique among the conservation integrals. The integral is 
divergence free for a material which admits (nonlinear) strain energy function. 
Furthermore /-integral has the same value for all open paths beginning on one face of the 
crack and ending on the opposite face (this assumes that there are no contribution to the 
integral from the crack faces, a condition which is usually met in most crack problems). 
This special path independent property under fairly general conditions has been 
advantageously exploited in the development of nonlinear fracture mechanics. For 
example, the above noted path independence allows a direct computation of the strength 
of crack tip singularities by evaluating the line integral in regions remote from the crack 
tip or along the remote boundaries. 

To the extent that the HRR singularity field (Hutchinson [8], Rice and Rosengren 
[9]) prevails over a distance which is larger than the fracture process zone and the zone of 
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the finite strain, the crack tip field can be said to be characterized by the value of J- 
integral and the onset of crack growth can be correlated with a critical value of J. The J 
value referred to here is evaluated on contour placed within a region dominated by HHR 
singularity. It is in this context that J is referred to as a crack tip parameter. In other 
words, the concept of a crack tip parameter is not tied to the property of global path 
independence. Nevertheless path independence has far reaching consequences, e.g. the 
value of crack tip parameter can be determined from the remote field. Within the 
deformation plasticity theory, J also has an energy based definition which has been used 
to derive useful formulae for its determination directly from the load-displacement record 
of a cracked body. In other words, while global path independence and the energy based 
definition are not essential to the concept of J as parameter characterizing the crack tip 
field, these features play important roles in an engineering fracture mechanics 
methodologies. 

The evaluation of crack tip contour integrals in numerical studies is a potential 
source of inaccuracy. To circumvent these numerical difficulties, Kishimoto, Aoki and 
Sakata [10-12] and Alturi, Nishioka, Brust and their coworkers [13-19] have restated 
crack tip integrals as path and area integrals. These so called path-independent integrals 
contain two terms: a path integral evaluated along a remote contour and an area integral 
evaluated within the area enclosed by remote contour. In three dimensional crack 
problems, these path and area integrals generalize to surface and volume integrals 
respectively. 

Accurate pointwise values of J along 3-D crack front have been obtained by Shih, 
Moran and Nakamura [20] and Nakamura, Shih and Freund [21]. Li, Shih and 
Needleman [22] have interpreted the method as an application of principal of virtual 
work and made the connection to the virtual crack extension technique of Parks [23, 24] 
and Hellen [25]. Under more restrictive conditions, a similar domain form was obtained 
bydeLorenzi [26]. 

A unified approach to the derivation of crack tip and associated finite domain 
integrals has been presented by Moran and Shih [27] using a general balance statement as 
the starting point. 
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One of the common methods used to evaluate the displacements and stresses 
appearing in the expression of /-integral is finite element method. 


1.3 OBJECTIVE OF THE PRESENT WORK 

The problem considered in the present work is given for study by the Computational 
Mechanics Section of the Reactor Safety Division (RSD) of Bhabha Atomic Research 
Centre (B.A.R.C.), Mumbai. The problem is concerned with the assessment of structural 
integrity of a typical, bimetallic welded pipe joint. The joint connects the to the pressure 
vessel of the reactor made of different material. So the accurate analysis of this joint is 
necessary. Numerical methods like FEM offers a viable and reliable solution technique to 
assess the severity of a crack developed in the welded pipe joints. 

The present work can be subdivided into two main parts, as detailed below, along 
with a brief description of each part: 

1. To develop an efficient Matlab code to evaluate the /-integral for three 
dimensional body having crack and validate it with known problems. 

2. To carry out a parametric study of a pipe joint having elliptical crack for different 
crack lengths and aspect ratios. 


1.4 STRUCTURE OF THESIS 

Mathematical formulation of J integral and its finite element implementation is described 
in Chapter 2. Finite element analysis package PATRAN/NASTRAN is used to generate 
the input data for finite element implementation of /-integral. Results and verification of 
the two dimensional and three dimensional problems are detailed in Chapter 3. Chapter 4 
deals with actual pipe problem. Results and discussion for the pipe problem are given in 
Chapter 4. Scope for future works has been discussed in Chapter 5. 
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CHAPTER 2 

MATHEMATICAL FORMULATION 


2.1 INTRODUCTION 

The study of a physical phenomenon usually involves two major tasks: Mathematical 
modeling of the physical process resulting in algebraic, differential or integral equations 
and solution of the resulting equation. 

While the derivation of the governing equations for most problems is not unduly 
difficult, their solution by exact methods of analysis is a difficult task. In such cases, 
approximate methods of analysis provide alternative means of finding solutions. 

The Finite Element (FE) method is one of the most popular methods of numerical 
analysis, in which, a given domain is represented as a collection of simple sub-domains, 
called finite elements, so that, it is possible to systematically construct the approximation 
functions over each element, needed in a variational or weighted-residual approximation 
of the solution of a problem. 

In this chapter, the mathematical formulation of crack tip integral, and FE 
formulations for a 3D static stress analysis have been presented. First we discuss the 
general derivation for crack tip integral. Subsequently 2D and 3D specialization will be 
presented. 


2.2 MOMENTUM BALANCE -CRACKTIP INTEGRAL 

2.2.1 General derivation of crack tip integral 

A general treatment of crack tip integral based on a variational form of momentum (or 
incremental momentum) balance is presented in this section. 


7 



Attention is restricted to small strains and the strain displacement relation takes 
the usual form 

£ ij = ( U u u j,i )/ 2 • ( 2J ) 

In the absence of body forces the equation of motion (or balance of linear 
moment) is written as 

( 2 - 2 ) 

where p is the mass density, a tj - <x jY is Cauchy stress and a superposed dot denotes the 
material time derivative. Taking the inner product of (2.2) with the velocity field, u t , and 
rearranging the resulting expression gives 

{a ji u i ) J =cXj iJ Ui+c7 Ji u iJ 

= pu t u t +aji itij (2.3) 

= f + W, 

where the stress work density, W, and kinetic energy density, T , at the material particle 
are given by 

w = | VijSij dt =>W = cry Sy = cry u tJ & 

T = | pu/Uj dt. =$T = piijU r 

Equation (2.3) is a differential form of mechanical energy balance and is valid for 
any mechanical response. 

The following fluxes are introduced for purpose of subsequent development 

<p j =cx ij u i , ij/ = (jy + t) (2.5) 

Here <pj is identified as the mechanical work flux vector and ys is the time rate of change 
mechanical energy density. The balance law (2.3) can be therefore written as 

<Pj,j = V- < 2 - 6 ) 
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The equation (2.6) is called general balance equation. Consider a small arbitrary 
volume V and a surface dV which surrounds the volume V . Integrating this expression 
over V and applying divergence theorem to the term on the left hand side yields 

\(Pjmj dS=\y/dV, (2.7) 

dV V 

where ntj is the outward unit normal to the surface dV . Here it is assumed that the fields 

are sufficiently smooth so that the divergence theorem may be applied. By applying 
divergence theorem to (2.3) the order of differential equation is reduced by one so the 
expression (2.7) can be viewed as a weak or variational form of linear moment balance. 
Since the variation field is velocity the expression (2.7) is also an integral form of 
mechanical energy balance. 

Let v,. be the instantaneous velocity of the surface dV . Using the Reynolds 
Transport theorem on the right hand side of (2.7) yields 

^(pjnijdS = — JV dV - ^y/Vjtrij dS. (2.8) 

3V dt y gy 

The result (2.8) is simply a representation of mechanical energy balance for the 
volume V . This result is now specialized to the crack propagation and for the purpose of 
clarity attention is focused on planer crack problems. 

Consider a two dimensional body with an extending crack oriented along the x x 
axis of a rectangular Cartesian coordinate system. The cracked body lies in the x x - x 2 
plane with the crack plane given by x 2 = 0 . The crack extends in its own plane along the 
x x axis with speed v . Define a vanishingly small crack tip region with a small contour T , 
which is fixed in the size and orientation with respect to the crack tip and is translating 
with crack tip at the speed v . As depicted in figure (2.1) the area bounded by the fixed 
material curve C 0 and the curve T is denoted by A(t) and is free of singularities. 
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Figure 2.1 Conventions at crack tip. Domain A is enclosed by F , C + ,C_, and C 0 .unit normal 
nij = rij on C + , C_ and C 0 , and ntj = —tij on T . 

Defining the contour C — C 0 +C + +T + C_ which encloses A(t) . Without loss of 
generality, the crack faces are taken to be traction free. Rewriting the left hand side of 
expression in (2.8) for two dimensional case 

Jtpj rrij dS = j ' (pj ntj dS + j (pj rrij dS - ^(pj rij dS + JV> 7 - m j dS 
dv c 0 c + r c. 

= J <pj rrij dS - j<P2 m 2 dS + j<Pj ntj dS + j<p 2 m i dS 
c 0 c + r c_ 

= jtpj v/ij dS + J ’<Pj ntj dS = j<Pj ntj dC + jVy ntj dT . 

c 0 r c 0 r 

Rewriting right hand side of the expression in (2.8) for two dimensional case 
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d 

d 


-\y/dV- vj mjdS = — fy/dV - j V v y . nij dS 

l y BV dt A C 


= -j^ ji// dV - Ji// Vj nij dS - JVv y . rrij dS 


= Yt \ ¥dV ~ dS 


= Tt$ y /dA -\y /v J m j dT - 


Substituting v y - = v 5 ^ and combining the left and right hand side as derived as two 
dimensional case of the expression (2.8) we get 


J \<Pj m j dC = ~ JysdA- f(<p. + if/ vS xj ) k ; dT. (2.9) 

c 0 at a r 


It is noted that the only non- vanishing velocity on the boundary of A(t) of consequence 
in (2.9) is Vj = v on T . The term on the left hand side of (2.9) is the rate at which the 
energy is being input into the body, the first term on the right hand side is the rate of 
increase of internal (mechanical) energy and consequently the last term is the 
instantaneous rate at which energy is being lost from the body due to flux through T . This 
quantity is denoted by F and is written for n y . = -m y . on F , 

F(rhS(p,+w-^ y )»jdr. ( 2 . 10 ) 

r 

There are two contributions to this flux integral. The first term represents the flux across 
T due to the material outside of T working on the material inside it. If T were a 
material curve this would be only contribution. The second term represents the 
contribution due to the flux of the material across T . The quantity G is defined as the 
limiting value of F(r)/v as the contour T is shrunk on to the crack tip, i.e. G is the 
energy released from the body per unit crack advance (per unit thickness). For this 
concept to have physical significance, the limiting value G must be independent of 
actual shape of T in the limit T ->• 0. In other words, the value of G must be path 



independent in the limit T — » 0 (i.e. in the crack tip region). We now consider the 
condition under which G is finite and path independent in the limiting sense described. 

Consider the closed path formed by two crack tip contours T, and T 2 and the 
crack face segments which connects the ends of the two contours. Application of 
divergence theorem to the closed contour integral yields 

f ( T 2 ) ~ ^( r i ) = \[(Pj + ¥ v ) y dA 

= \(<Pjj + y / j vS \j) dA ( 2 - n ) 

A\ 2 

= + v,\ v ) dA 

A\2 

where A n is the area enclosed by the closed contour. By using (2.6) in (2.1 1) we obtain 

^(r 2 )-F(r,)= jV + vvo)<£4. (212) 

A\2 

If the integrand in (2.12) is o(l/r 2 ) where r = *Jx x 2 + x 2 2 then J p(r 2 )-F(r,) = 0 as 

r,,r 2 — > 0 and path independence of F in the crack tip region (r — > 0 + ) is established. 
It is noted that this condition is satisfied if any field quantity / satisfied the following 
relation. 

f + v /,i = o(f tX ) as r -> 0 + , (2.13) 

where r is the radial distance from the moving crack tip. The condition (2.13) can be 
interpreted as a condition for locally steady state behavior. In particular if, y/ satisfies 

(2.13) then y/ is of the order of 1/r, we have 

iff + vy/ x =o(l/r 2 ) asr-»0 + , (2.14) 
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and the path independence in the crack tip region is assured, it follows that on choosing 
T to be circular contour, (2.10) yields a finite value of F if the integrand of (2.10) is of 
order (l /r), i.e. 


(<Pj +ifA>8 X j)t* Aj(0)/r + o(l/r) as r— »0 + , (2.15) 

for some Aj(o) condition (2.14) may also be viewed as an integrability condition. 

For field variables, <p andi//, satisfying a balance law of (2.6) and for which local 
steady state conditions (2.13) and (2.15) hold, the integral given by (2.10) is path 
independent within the local crack tip region. It is noted that (/ + yf, ) is the time 
derivative of / with respect to a coordinate system which translate with the propagating 
crack and zero for a steady state condition. Thus the right hand side of (2.14) vanishes 
identically under the steady state crack propagation conditions and (2.10) provide a 
globally path independent integral for this class of problem for any material response. 

With the understanding that the limiting value of fluxF(r), denoted by F is 
independent of shape of contour F as it shrunk on to the crack tip. The relation for F is 
written as 

F = lim |(^ y - + v y/S XJ )« y - dT. (2.1 6) 

r_>0 r 

A general representation for energy release rate can be derived using the general flux 
tensor H^- . With its *, component denoted by H ly , the general representation of energy 

release rate can be written as 


J = — = lim [h,.- rijdr. (2.17) 

v r->oJ J J 

It follows from (2.16) that H, ; is defined by asymptotic relation 

Hiy ^ (<Pj + vyS xj )/v 

i a 


as r-¥ 0 + . 


(2.18) 



The above relation is not written as an equality since the local steady state condition 
(2.13) will necessarily be invoked in defining H ly - . Assuming that the flux field in H ly - 

are sufficiently smooth and noting (2.6) and the asymptotic conditions (2.14) and (2.18), 
the condition for J to be finite can be restated as 

H ly - « Aj {d)/ r + o(l/ r) and 6, = H ly . y - =o(l/r 2 ) asr-»0 + . (2.19) 

It is noted that with (2.6) as the starting point no additional physical principle has been 
invoked to obtain result embodied in (2.16) and (2.17). 

2.2.2 Mechanical energy integral 

When cpj and if/ are identified by (2.5) the flux integral (2.16) has the explicit form 

F = lim j[(r + T)vS,j +a„ dT. (2.21) 

~T 

The above relation is well known integral expression for crack tip mechanical energy 
flux. To obtain the corresponding result for energy release per unit crack advance J 
commonly referred to as energy release rate, we make the substitution u i &-vu iX (i.e. 

local steady state condition (2.13) has been applied to the velocity field) and rearrange 
(2.20) to get 

y=lta \[(W + T)S XJ -^u^ nj dr, (2.22) 

u r 

by formulation (and necessity) the loop F translates with the moving crack tip. Under the 
transient conditions or general material response, the energy release rate is given by the 
crack tip limit in (2.21), although the shape of the contour as shrunk on to crack tip is 
arbitrary. 

Comparing (2.21) with expression (2.17), the integrand may be written more 
compactly as X ] component of general energy tensor, 

H, ; =(H' + T)S, l -a l u u , 


(2.23a) 



the general energy tensor can be given by 

H 9 =(W + T)S 9 -0 g u u , (2.23b) 

2.3 A DOMAIN INTEGRAL EXPRESSION FOR THE ENERGY 
RELEASE RATE 

2.3.1 Two-dimensional specialization 

Consider a two-dimensional body and direct attention to a line a line crack along 
the x x axis. For arbitrary crack tip motion (but restricted to crack advance along the x x 

axis) under the dynamic conditions the energy release per unit crack advance is given by 
equation (2.22) 

J=\im\iw + T)S lJ -a IJ u, 1 \n l dC, 

r 

where W and T are the stress work and kinetic energy densities respectively, cr ( y and w ( . 
are the Cartesian components of the stress and displacement, n i is the unit vector normal 
to T and dC is the arc length as depicted in figure (2.2). 

In the context of linear elastic material response, the result (2.22) is proposed by 
Atkinson and Eshelby [28] and independently derived by Kostrov and Nikitin [29] and 
Freund [30]. As noted by Nakamura et al. [31], “the general result (2.22) underlies 
virtually all crack tip energy integrals that have been defined and applied in fracture 
mechanics in the sense that each may be derived from (2.22) by invoking the appropriate 
restriction on material response (through W ) and on the crack tip motion”. It is also 
noted (2.22) was obtained from purely mechanical consideration. With usual 
reinterpretation of the field quantities with reference to undeformed configuration, (2.22) 
also hold in the context of finite deformation. 
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Figure 2.2 Simply connected region A } enclosed by contour C = C 2 + C 3 — C x + C 4 


Under quasi-static condition, T = 0 , and (2.22) reduces to 

J=pmffc'S,j-<T j/ uJn J dC. (2.24) 

r 

For deformation theory solid, and in the absence of body forces, thermal strains and crack 
face traction the integral in (2.24) is path independent and its value does not depend on a 
limiting process in which the crack tip contour F shrunk onto the crack tip. The above 
result is due to Rice [1], For this important case, the contour integral value for (2.24) can 
be accurately and readily extracted from the finite element field remote from the crack 
tip. It is noted that Eshelby introduced the conservation integral (2.24) in the context of a 
force on a dislocation or a point defect; under the isothermal conditions the bracketed 
quantity is x, component of energy momentum tensor H (see 2.23a). 

The equation (2.24) can be rewritten as [20, 22] 
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J = [[^ij u ij - w5 \ iK <h dC ~ [ 3+Ci 


cr 2J u jl m 2 q ] dC. 


(2.25) 


Here the closed curve C = C 2 + C 3 - C, + C 4 , is the sufficiently smooth function in 
the region enclosed by C that is unity at C, and vanishes on C 2 and ntj is outward 
normal as depicted in Fig. 2.2; m 1 = 0, m 2 = ±1 on the crack faces, and trij = -rij on 
Cj . The last term in (2.25) vanishes for traction free surfaces. 


Within the framework of small deformation theory, the strain energy density can be 
written as the 




(2.26) 


Applying divergence theorem to closed contour integral in (2.25), we get 

J= [, u j-' ~ w s " L dA - U a2j Uja mi qx dC 


(2.27) 


where A, is the simply connected (area) domain enclosed by C. Now invoke the 
equilibrium condition and differentiating (2.26) with respect to x, , we get 


<r 9 j+fi=0> „ 


(2.28) 


to obtain the desired domain expression for the energy release rate 


J= L, “ W s \ i fci ./ ~ ft “u <h} dA ~ h «/,i dC. (2.29) 


Here is the body force per unit volume and t. is the traction on the crack faces. In the 
absence of body force and crack face traction, (2.29) reduces to 

■'”lir,Uj ) -WS u ]q u dA, (2.30) 

JA ! 

which is energy release rate expression given in [22,26]. Equation (2.30) is domain 
independent in the sense that any annular (area) domain can be chosen for the purpose of 
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evaluating / . In this case the choice of domain is often dictated by convenience and mesh 
design. With the presence of thermal strain, body force and/or crack face traction, the 
domain integration for (2.30) must include the crack tip region (r -» 0 + ) . 

The function q x can be interpreted as imposing on the material points on C, a unit 
translation in x, direction while requiring the material points on C 2 to remain fixed. The 
material points within A, are displaced according to any smooth interpolating function. 
In the context of a finite element model, C, is the arc formed by crack tip node/nodes. By 
letting C, coincide with the boundaries of different rings of the elements surrounding the 
crack tip, J can be evaluated on the alternate domains. 

2.3.2 Three dimensional formulation 

We consider a three dimensional crack front with a continuously turning tangent. Imagine 
that the plane shown in the Fig.2.2 is now the x, - x 2 plane of the local orthogonal 

Cartesian system depicted in Fig. 2.3a. Asymptotically, as r — » 0 + , plain strain condition 
prevail so that the three dimensional fields approach the (plain strain) two-dimensional 
fields at the crack front. Thus (2.24) defines the pointwise energy release rate which can 
be denoted by J(s) 


1R 



(PERPENDICULAR 


*2 



Figure 2.3 (a) Definition of the local orthogonal Cartesian coordinates at the point s on the crack front; 
the crack plane is in the X, — X 3 plane, (b) The function 8 l(s) can be pictured as a virtual crack advance 
in the direction normal to the crack front in the crack front. 


Let Sl(s) denote the local crack front advance at the point s in the direction normal to 
the crack front and in the plane of the crack, and let ds denote the elemental arc length 
along the crack front as depicted in the Fig.2.3(b). Recalling the definition /-integral, 
which states that it is the energy release rate per unit crack length advance i.e. 
/ = -dll/da where -dll is decrease in potential energy at a point on crack front. 
Then to within first order in 8 1 

-8U = JAa= | J(s)8l(s) ds. (2.31) 

where L c denotes the line segment of crack front under consideration and 81 is an 
arbitrary crack front advance. Here -£IT is the decrease in potential energy and 7 is 
the energy released per unit of the finite segment of crack advance. 

By separately advancing various segment of the crack front, the arc weighed 
pointwise J(s) can be computed along a three dimensional crack front. 

To fix ideas and minimize nonessential manipulations we focus our attention on a 
notch with notch thickness h as shown in Fig. 2.4 and argue (heuristically) that h 0 is 
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the sharp crack configuration of interest. The surface of the notch consists of faces S A 
and > S B , with normals along x 2 , and a face with a normal in the x } - x 3 plane. Now let 
the notch face S t with normal m k in the x, -x 3 plane, advances A al k in the x k 
direction, i.e. 

Aa.l . m = 8l{s) or A al k m k =Sl{s) (2.32) 

Furthermore restrict l k to be a function of and *3 only; thus / 2 s 0 . Using (2.23b) and 

(2.32) in (2.31), we obtain 

J Aa = Aa [[ a ij u j,k- ws k\h m i ds - -^ 2 - 33 ) 

When A al k corresponds to a translation, l k can be taken outside of the integral sign and 

(2.33) gives J=l k J k where under isothermal conditions J k are the translational 
conservation integrals (Budiansky and Rice [7]). 
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Figure 2.4(a) Schematic of body and volume V in the Jtj - x 2 plane containing a notch of thickness 
h ;note that m = — n on S t . (b) Schematic of notch face when the function A al k is interpreted as a 
virtual advance of a notch face segment in the direction normal to x 2 

Consider the simply connected volume, V, enclosed by surfaces S n S + ,S~ and^ as 
depicted in Fig. 2.4(b). We introduce the function q k defined by 



< h = 


(2.34) 


[l k on S. 

[0 on S' 

together with the requirement that q k be sufficiently smooth in the volume V . Using 
(2.34), (2.33) can be rewritten in the form of an integral over the closed surface 
S (. S = iS, +S + +S~ -S t ) plus the crack face terms, 

j i [°v u j,k - W 5 ki K q k ds - a 2j U M m 2 q k ds. (2.35) 

To arrive at (2.35), following results have been used m x = 0, m 3 = 0 (m 2 = ±1), on the 
crack faces; it may be noted that l 2 =q 2 =0 everywhere. 

Application of divergence theorem to the closed surface integral in (2.35), and use of 
(2.28) yields 

J = I { [°> u j,k ~ W S b - ft u u <lk}dV- t, u ik q k ds. (2.36) 

Now let h — > 0 to obtain the desired to obtain desired expression for the energy decrease 
when a local segment of crack front advances by A al k in its plane. In the absence of 
body force and crack face traction (2.36) reduces to 

(2-37) 

The domain expressions (2.36-2.37) give the energy release by the body per unit crack 
advance over a finite segment of the crack front. To a first approximation the pointwise 
energy release rate J(s) may be assumed to constant over the small length L c of crack 
front and can be brought out of the integral sign in (2.31) to yield, 

A») = , (2.38) 

A a ^ l k ( 5 ) m k (j) ds 

where the term in the denominator is the increase of crack area due to the virtual crack 
advance at point 5 on the crack front and in Jis evaluated over the volume which 
contains the length L c . 
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2.4 FINITE ELEMENT FORMAULATION FOR THE VOLUME 
INTEGRAL METHOD 

2.4.1 Two-dimensional implementation 

For isoparametric elements, the coordinates (jc,, x 2 ) in the physical space and the 

displacements («,, u 2 ) for 8 noded element are written as 
* 8 

x i = 2X X Kt ii, = y Z N xU iK , i = 1, 2 (2.39) 

K=l K - 1 

where N K are the bilinear shape functions, X ik are the nodal coordinates and U, K are the 
nodal displacements [32]. 

Consistent with the isoparametric element formulation, q x can be interpolated within an 
element as 


<?> =IXG„ (2.40) 

7=1 

where Q x , are nodal values for the I th node. From the definition of q x , if the I th node is 
on C, , Q x , = 1 , whereas if the I th node is on C 2 , Q u = 0 . In the area between C, and C 2 , 
Q XI can be taken to vary between 1 andO . Using (2.39), (2.40) and the chain rule, the 
spatial gradient of q x within an element is given by 


dffi _ v-iy- dq 
dxj %£(d Vk dI^ u 


(2.41) 


where drjjdxj is the inverse Jacobian matrix of the transformation (2.39). 


With 2x2 Gaussian integration, the discretized from the domain expression for the 
energy release rate (2.29) for plane strain and plane stress problem is 
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j- Z I 

all p - 1 
elements 
in A | 


du 


dx x 


L -WS „ 


dgi f du, 

dx, Ji dx x qx 


det 


' a *.' 1 


3 77 




\ w • k J 


- 2 

all edges 
onC 3+C4 


t. 


3«; 

3x 


■q } Vw. 


(2.42) 


Here the quantities within { are evaluated at the 4 Gauss points and and w are 

the respective weights. The crack face traction contribution can be evaluated using the 
equivalent nodal forces and spatial gradients of the nodal displacements or using standard 
integration procedure for boundary pressures [32]. (Note that the integration is carried.out 
only on element edges on C 3 and C 4 with non-zero traction). 

2.4.2 Three dimensional implementation 


Consider a segment of the crack front, e.g., as shown in Fig 2.4. The limit h-+0 has 
been taken to obtain the crack line identified by L c . In the subsequent discussion L c will 
be taken as the line connecting the nodes M -1, M and M + 1 as shown in Fig. 2.5, and 
identify the volume V with the collection of elements which contains the lineZ, c . 

Consistent with the isoparametric formulation q i within an element can be interpolated 


as 


*-1.2.3, (2-43) 

/=! 

where N, are the trilinear shape function and Q u are the nodal values for the I th node. 
From the definition of q k (see (2.34)), Q u = 0 if the I th node is on 5, . For nodes inside 
V , Qu is given by interpolating between the nodal values on L c and S x ; furthermore Q 2! 
vanishes identically. 
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Figure 2.5 Schematic of the finite element mesh in the x, - jc 3 plane illustrating the basis function for an 

“interior” node and “end” node along the crack front 

Following the standard manipulations, 


dffi _y 'yi dMy g?y 
dXj khdv k ~d^~ in 


(2.44) 


where rj k are the coordinates in the isoparametric space. Employing 2x2x2 Gaussian 
integration, the discretized form of the (2.36) is 


£ £- 


all p=l 
elements 
in V 




QUj 


v ax* 




S.X; 


‘"/I 


du i 

dx t 


(Ik 



( dx t ^ 

det 

k 


{ 


w„ 


- £ • 
all edges 
011CJ+C4 


t; 


du. 


' dx. 


■ <lk- 


(2.45) 


where J is the energy released for a unit virtual advance of a finite crack front segment. 
Here the quantities within { } p are evaluated at the 8 Gauss points and w p and w are 
the respective weights. As discussed previously, the crack face traction contribution can 
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be evaluated using the equivalent nodal forces and the spatial gradients of the nodal 
displacements, or using the standard integration procedure for boundary pressure [32]. 

The relationship between J(s) and J can be established by substituting (2.32) into 
(2.31) to get 

| J(s)lj(s)mj(s) ds = J. (2.46) 


Now introduce piecewise continuous function for J(s) and lj(s)mj(s). 

Let J K denote the value of J(s) at the K‘ h node on the crack front and ¥*($) be 
piecewise linear functions defined along the crack front. Then 

;(.)=!;«(>). < 2 - 47 > 

K = 1 


where N is the total number of nodes on the crack front and 


Jl at the K th node 

[O at all other nodes 


The nodal value J K can be obtained using (2.38) 


r [¥*(*)*’ 
*0 


(2.48) 


(2.49) 


2.5FE FORMULATION FOR 3D STATIC ELASTIC STRESS 

ANALYSIS 

Consider an arbitrary three-dimensional solid body Q with volume V and boundary 
r as shown in Figure 2.6. 
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r 


Figure 2.6 Domain and Boundary Conditions for Stress Analysis. 

With reference to the Figure, T u represents the boundary on which displacement 
is specified. Traction t is specified on the boundary F, and P c represents the 
concentrated forces acting on the body. 

The principle of virtual work states that for a virtual perturbation in the 
displacement field at equilibrium, the virtual change in the internal strain energy (SU) 
must be balanced by an identical change in the external work (8W ) due to the applied 
loads, that is: 

SU = SW , (2.50) 

or. (2-51) 

e e 

where, SU C and SW e are the contributions of the finite element e to the virtual change in 

strain energy and external work respectively. 

Now, 

SU e = \{Ss} T {cr}dV e , (2.52) 

where, V e represents the volume of the element e , 

{e } = Strain vector = [s„ s w s lz , 

{<r}= Stress vector = [cr^ <7 a cr^ crj^ , 
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Assuming an isotropic linear elastic material behaviour, (2.52) may be expanded as: 


eu. = \{ss} T [D\{s\dv,, 

v e 

where, \p\ is the material property matrix given by: 


(2.53) 


M- 


(l + uXl-2u) 


l-o 

V 

V 

0 

0 

0 


0 

1 -v 

V 

0 

0 

0 


v 

V 

1-0 

0 

0 

0 


0 

0 

0 

l-2o 

2 

0 

0 


0 

0 

0 

0 

l-2o 

2 

0 


0 

0 

0 

0 

0 

l-2o 


where, E is the Young’s modulus of the material and o is the Poisson’s ratio; 
The strains are related to the displacements as: 


(4= 


du dv dw 
dx dy dz 


( du chA ( dv ( du dw'' \ 


■+ — 
dy dx 


[dz + dy y 


dz dx 


(2.54) 


where, u , v and w represent the displacements of any arbitrary point within the element, 
in x , y and z directions respectively, and are interpolated from the nodal displacements 
through the element shape functions as: 


r \ 

U 

V 

w 


=Mki. 


(2.55) 


where, [iV] is the element shape function matrix and {u e } is the element nodal 

displacement vector. 

From (2.54) and (2.55), 

< 2 - 56 ) 

where, [5] is the strain-displacement matrix based on the derivatives of the element 
shape functions. 
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Combining (2.53) and (2.56) and noting that \u e ) does not vary over the volume of the 
element: 

XJ. - {*#.}' PY[D][B]dV,{U'}dV , . (2.57) 

V e 

Now, 

W. = {Sh.Y + {Su,Y{p:}, (2.58) 

f 

where, T, e is the area of the element that falls on T, , 

}= Traction vector on T, = [/* t e y t\ ] r , 

{P c e } = Concentrated nodal force vector. .. . 

Finally, substituting (2.57) and (2.58) into (2.51) and writing the Equation on element 
basis to give: 

iK J. [Bj[D][B\dV, {u, } dV, f = %({Su. Y | M r If }< + {Su, f {p ; }) (2.59) 
Now, since {S u e } r is arbitrary, (2.59) reduces to: 

EMW-ZkblM+Ifc*}. < 2 - 60 > 

e e e e 

where, [F e ]= J[#] r [l)][5]rfF e = Element stiffness matrix, 

v. 

}= JpV ] r {f e = Element traction load vector, 

r, e 

JWM* m }^K = Element thermal load vector. 

K 

The element matrices and load vectors of all the elements in the domain are assembled 
together to get the global stiffness matrix [i^] and the global load vector {f} such that: 

where, \u g } is the global nodal displacement vector. 
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CHAPTER 3 

VERIFICATION AND VALIDATION OF RESULTS 


3.1 INTRODUCTION 

The closed form solutions for the stress intensity factor (SEF in mode-I) of simple test 
specimens like Compact Tension (CT), Single Edge Notch Bend (SENB) specimen _and 
Centre-Cracked-Plate (CP) are known. In elastic region, a relation between J and K, 

exists. So the value of J can be evaluated accurately for a simple elastic problems. 
Thereby, the code can be verified with the help of these problems. First, the code for 2D 
case is verified for the SENB, CT and CP specimens under plane strain. Subsequently the 
code for 3D case is verified for solving the same set of problems with thickness is taken 
into consideration. 


3.2 DESCRIPTION OF ANALYSIS 

The stress analysis is performed using the PATRAN/NASTRAN FEM package. To 
determine the stresses at the Gauss points, first the nodal coordinate and displacement 
data is extracted from the analysis and then, with the help of these displacements, stress 
tensor, strain energy and finally J integral are evaluated using a Matlab code. 

For 2-dimensional problems, 8-noded isoparametric quadratic elements have been 
used and the problems are analyzed in plain strain conditions. For three dimensional 
cases, 8-noded trilinear brick elements have been used. 



3.3 VERIFICATION OF CODE FOR TWO DIMENSIONAL 
PROBLEMS 

The results obtained from the Matlab code are verified for three different problems: 
SENB, CT and CP specimen. For two dimensional cases, the domain for ./-integral is an 
area bounded by two curves around the crack as shown in Fig. 2.2. As discussed 
previously, in the absence of crack face traction and thermal strain, the inner boundary of 
the domain need not be too close to the crack tip. In this way, the stress field near the 
crack tip (which may contain some inaccuracies) does contribute to the calculation of J 
integral. In fact, the variation in J values for two dimensional cases is less then 0.1% for 
the domains well within the specimen (far from the specimen boundary). 

3.3.1 SENB Specimen 

Geometry of the specimen is shown in Fig. 3.1. The specimen is symmetric with respect 
to the plane in which the notch lies. We analyze only half of the specimen taking benefit 
of the symmetry. The mesh with boundary conditions is shown in Fig. 3.2 in which 1 
represent constrain in X direction and 2 represent constraint in Y direction. 



T 

w 

i 
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8 fixed inner 
boundary 


Crack tip 


Figure 3.3 Domain selection. 
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The analytical expression for the Stress Intensity Factor K, for the SENB 
specimen is given by [34] 


K, 


PS_ 

BW V2 


/(*); 



(3.1) 


where, B = Plate thickness 


and 


m= 


3a l/2 [l.99-a(l-q) (2.15 -3.93a + 2.7cr 2 )] 
2{l + 2a)(l-af 2 


(3.2) 


In the elastic region the relationship between /-integral and K, for the plane strain case 
can be written as 


/ = 



(3.3) 


where E and v are elasticity constants. 

The model dimensions and loading for analysis are taken from [33], and are the 


following. 

Span 5 = 4 W = 101.6 mm 
Width W = 25.4 mm 
Crack Length, a = 12.7 mm 

Plate Thickness B = 25.4 mm (3.4) 

a/W = 0.5 

Applied load P = 802.5 N 
Analysis load = P/2 = 401.25 N 


To demonstrate the domain independence of/ -integral, it is evaluated on a number 
of domains. A set of domain can be created as shown in Fig. 3.3. Domains are obtained 
by adding layer of elements to a smaller domain. Several contours are taken in each 
domain (see Table 3.1) by making the inner boundary layer (pushing it closer to the outer 
boundary). We evaluated the /-integral corresponding to each of these domains and 
results are tabulated in the Table 3.1. This table contains values of the /-integral for 9 
sets of domains shown in Fig. 3.3. In the first set there are 14 contours and in the 
subsequent sets the number of contours decrease. This is because in each set the outer 
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boundary of the outermost domain (domain in 14 th row) is fixed and the inner boundary is 
moving outwards so each time the number of domain is decreased by one. 


Table 3.1 Values of J -integral (in N / m ) for different sets of domains corresponding to different inner 
boundaries for a theoretical 7-integral value of 1 8 .40 N /m ■ 


hi 

1 

2 

3 

4 

5 

6 

7 

8 

9 


18.7284 










18.7307 

18.7318 









18.7334 

18.7348 

18.7372 








18.7347 

18.7356 

18.7370 

18.7367 







18.7356 

18.7363 

18.7373 

18.7373 

18.7378 






18.7360 

18.7364 

18.7371 

18.7371 

18.7372 

18.7367 




7 

18.7362 

18.7364 

18.7369 

18.7368 

18.7368 

18.7364 

18.7361 



8 

18.7364 

18.7365 

18.7369 

18.7368 

18.7368 

18.7365 

18.7364 

18.7366 


D 

18.7365 

18.7365 

18.7368 

18.7368 

18.7368 

18.7365 

18.7365 

18.7366 

T 8.7365 

19 

18.7365 

18.7365 

18.7368 

18.7367 

18.7367 

18.7366 

18.7365 

18.7366 

18.7372 

n 

18.7365 

18.7365 

18.7367 

18.7367 

18.7367 

18.7366 

18.7365 

18.7366 

18.7371 

irf 

18.7365 

18.7366 

18.7368 

18.7368 

18.7368 

18.7367 

18.7367 

18.7367 

18.7375 


18.7398 

18.7399 

18.7403 

18.7405 

18.7407 

18.7409 

18.7413 

18.7419 

18.7431 


18.7845 

18.7863 

18.7886 

18.7910 

18.7941 

18.7977 

18.8025 

18.8088 

18.8173 


From the table, it is clear that the difference in the values of the J -integral for 
different domains is very less (less than 0.1%) for the domains well within the specimen. 
This error can completely attributed to numerical rounding off during evaluation of 
integration and loss of accuracy in data transfer from PATRAN. As we move downward 
in the table, the outer boundaries of the domains get closer to the specimen. Some portion 
of the outer boundaries of the domains corresponding to 14 th row coincides with the 
specimen boundaries where the stress values are not accurate, and hence the deviation in 
J value in these cases from the mean value of 1 8.7364 N/m is more. 

From Equation (3.1), (3.2) and (3.3) we can evaluate the value of J -integral in the 
elastic region for SENB specimen with dimensions and loading given by (3.4). The value 
of J -integral is found to be 18.40 N/m . The same specimen has been analyzed by [34] 
and the value of J is found to bel8.46 N/m . The error in evaluating 7-integral from the 
code is 1.82 % (with respect to /obtained from the formula). 
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3.3.2 CT Specimen 

The geometry of the CT specimen is depicted in Fig. 3.4. The specimen is being 
symmetric with respect to the plane in which the notch lies. We analyze only half of the 
specimen using the symmetry. The mesh with boundary conditions is shown in Fig. 3.5. 



Figure 3.4 Compact tension (CT) specimen. 





Figure 3.5 Mesh (with B.C. and loading) for CT specimen. 


The analytical formula for Stress Intensity Factor K, for CT specimen is given by 


[34] 


K, 


P 

BW' 12 


/(«); 



(3.5) 


where, B = Plate thickness 

„ N (2 + a)[o.886 + 4.64 a -13.32a 2 +14.72 a 3 -5.6a 4 ) 

/(a) “ (1 -a) 312 

Here E and v are elasticity constants. 

The model dimensions and loading for analysis are taken from [33], and are: 


(3.6) 


Crack length a = 14.9 mm 
Width W = 25 mm 

Thikness B = 25 .2 mm (3.7) 

a/W = 0.596 

Applied load P = 3806 N. 


We evaluated J -integral on the same sets of domains and found the similar trend for this 
specimen too. The results for few domains are tabulated below. 


Table 3.2 Values of J - integral evaluated on different domains for CT specimen. 


Domain 

J -integral (in N/m ) 

1 

710.8870 

2 

710.9680 

3 

710.9682 

4 

710.9336 

5 

710.9434 

6 

710.9434 

7 

710.9494 

8 

710.8828 


Value of J -integral as obtained using (3.5), (3.6) and (3.3) using the dimensions and 
loading given by (3.7) is 723.2 N/m. The error in evaluating /-integral from code is 
1.7%. 


3.3.3 Centre-Cracked Plate under uniform tension 

The geometry of the center-cracked plate is depicted in Fig 3.6 again, using the symmetry 
of the problem, we analyze only one fourth of the model. The mesh with boundary 
condition is shown in Fig.3.7 in which 1 represent constrain in X direction and 2 
represent constraint in Y direction. 

The analytical expression for Stress Intensity Factor K, for CP specimen is given 
by [35] 

K, = <xjn~a f(cc); a = -^r, ( 3 - 8 ) 

W 

and /(a) = 1.0 + 0.128 a -0.288 a 1 +1.523 a 3 (3.9) 
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or 


Figure 3.6 Centre-Cracked Plate under uniform tension 
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Figure 3.7 Mesh for Centre-Cracked Plate 
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The model dimensions and loading for analysis are as follows 
Crack length a - 0.25 m 
Width W = 0.5 m 
a/W = 0.5 

Applied traction t = 80 MPa 


(3.10) 


We evaluated /-integral on the same sets of domains and again found the similar 
trend for this specimen too. The results for few domains are tabulated below. 


Table 3.3 Values of / - integral evaluated on different domains for CP specimen. 


Domain 

/ -integral (in N/m) 

1 

30633.0 

2 

30633.0 

3 

30632.6 

4 

30632.2 

5 

30632.6 

6 

30632.4 

7 

30632.6 

8 

30632.6 


Value of /-integral as obtained using (3.8), (3.9) and (3.3) using the dimensions and 
loading given by (3.10) is 30451 N/m . The error in evaluating /-integral is 0.6 %. 
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3.4 VERIFICATION OF CODE FOR THREE DIMENSIONAL 
PROBLEMS 

The 3-D Matlab code is used to obtain the /-integral values in the same specimen used in 
sec 3.3. For three dimensional case the domain is a volume bounded by the surfaces 
which surround the crack-front, such that the inner boundary of the domain coincides 
with the crack front. In this way the stress field near the crack front (which contains some 
inaccuracies) induces some error into the value of /-integral is not as accurate as obtained 
for 2D case. In fact, the variation in / values for first domain, which is nearest to the 
crack front (and consists of layer of single element around crack front), is highest. As the 
number of elements in the domain increases, the value of /-integral approaches to the 
theoretical value. 

3.4.1 SENB Specimen 

The geometry and the mesh are shown in Fig. 3.8(a) and 3.8(b) respectively. We have 
taken the same specimen and loading as in case two dimensional case (see sec. 3.2). If we 
consider the thickness of the specimen in z direction then specimen is symmetric with 
respect to the x-y plane. The specimen is also symmetric with respect to the plane in 
which the crack lies. Hence we analyze only one fourth of the specimen. 
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(a). 



(b) 

Figure 3.8 (a) 3-d SENB Specimen(b) Mesh For SENB Specimen. 

Fig. 3.9 shows the variation of ./-integral (in N/m ) along the crack front. The pointwise J- 
integral is evaluated over 10 domains. 
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normalized distance 


free surface 


Figure 3.9 Variation of point- wise 7-integral along the crack front for SENB specimen. 


From the Fig. 3.9, it is clear that in the middle of the specimen, where there is 
plain strain condition prevails, the value of /-integral approaches a value slightly above 
(about 10%) that of plain strain case (18.7 N/m). As we move towards the free surface 
the value of /-integral decreases. Similar trend has been obtained for SENB specimen by 
Shih et. al.[20]. For the first domain, the deviation of pointwise / value is largest in the 
plane strain region in the middle of the domain. As we increase the elements in the 
domain by adding more layers of elements, the pointwise /-integral value approaches to 
plain strain value in the middle of the domain. From the Fig. 3.9 it is also clear that the 
domain with only three layers of element is sufficient to get approximate value of /. 


3.4.2 CT Specimen 

We have taken the same specimen and loading as in two dimensional case. Again only 
one fourth of the specimen is analyzed on account of the symmetry Fig. 3.11 shows the 
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is evaluated o\er 10 domains. The mesh of CT specimen used in analysis is shown 
Fig. 3.10. 



Figure 3.10 Mesh for 3-d CT specimen 



Figure 3.1 1 Variation of pointwise 7-integral along the crack front for CT specimen 

Similar trend as in case of SENB specimen is obtained here. In the middle of the 
specimen, the value of 7 close the plain strain of llON/m and as we move toward the 
free surface the value of ./-integral decreases. The value of ./-integral again is slightly 
greater (about 1 0%) than the plane strain case (71 0 N/m ) in the middle of the specimen. 


3.4.3 Centre-Cracked-Plate under uniform tension 

Fig. 3.10 shows the variation of point- wise 7-integral (in N/m) along the crack front. The 
point-wise 7-integral is evaluated over 10 domains. 
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Figure 3.12 Variation of pointwise 7-integral along the crack front for Centre-Cracked-Plate 

The trend obtained for this case (Fig. 3.12) is slightly different than the trends 
which have been obtained previously. As we move in thickness direction from the middle 
thickness of the plate to free surface the value of point-wise 7-integral increases and then 
decreases near the free surface. The value of 7-integral approaches a value slightly (less 
than 5% for higher domains) above the plane strain value (30630 N/m ) in the middle of 
the specimen. 
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CHAPTER 4 

RESULTS AND DISCUSSIONS 


4.1 INTRODUCTION 

Welded joints are used to join a pipe to the nozzle of the pressure vessel. Because of the 
welded joint at the intersection, there always exists a possibility of an initial flaw. Due to 
the high stresses and corrosion, an existing flaw may grow into a surface crack. Though 
the actual shapes of such cracks are not simple, it has been observed that they are usually 
almost semi-elliptical. Hence, for the purpose of numerical analysis, such cracks can be 
approximated by semi-elliptical profiles. In this Chapter, the effect of the tensile and 
bending loads on such semi-elliptical cracks in pipe is presented. 


4.2 FINITE ELEMENT ANALYSIS OF CRACKED PIPE 
4.2.1 Solid model and mesh for analysis 

The solid model of the pipe is shown in Fig. 4.1. Figures 4.2 and 4.3 show the 
X - Y cross section of pipe, in the middle of the length inZ direction. Crack is assumed 
to have symmetrical shape with respect to the plane 1 and it lies in the “crack plane” (in 
the middle of the pipe) as shown in Figs. 4.2 and 4.3. Taking the advantage of this 
symmetry we analyze only half of the pipe. The mesh for the pipe has been generated by 
a FORTRAN code. One of the meshes generated through the code has been shown in 
Figs. 4.4 and 4.5. A closer view of the mesh is shown in Fig. 4.6 and the crack-front is 
shown in Fig. 4.7. The mesh is imported into PATRAN/NASTRAN and elastic analyses 
are carried out. For the finite element analysis, 8 noded brick (hexahedral) elements have 




Figure4.1 Isometric view of the solid model of pipe. 



Figure 4.2 Cross-sectional view of pipe with crack. 
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Figure 



4.3 Closer view 





Figure 


4.4 Isometric view 


of a typical mesh- 
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Figure 4*6 Closer view of mesh near crack front 


crack 

front 



Figure 4.1 Mesh near crack- front. 


4.2.1 Geometric and material properties 

Geometric data of the pipe provided by BARC Mumbai is as follows: 

• Nominal Diameter of the pipe: 261mm 

• Thickness of pipe: 32 mm 

Pipe has a semi-elliptical crack which has semi major axis a and semi minor axis b as 
shown in Fig 4.3. The semi major axis is along the X direction. A pipe of 50 mm 
diameter is also analyzed for the purpose of determining the effect of pipe diameter on 
the results. Material data of the pipe for elastic analysis are taken as follows: 

• Young modulus: 210 GPa and 

• Poisson’s ratio: 0.30 

4.2.2 Loads and Boundary conditions 

Since the pressure does not have any effect on opening of the crack, model is analyzed 
for only two different loadings: tension and bending. The loads are as follows. 

Tension: 80 MPa 


50 


Bending Moment: 25 kNm 


The problem being elastic, the exact value of the load is not important. Figure 4.1 
shows different planes of the solid model on which the load and boundary conditions are 
applied. As discussed earlier, the model of pipe is symmetric with respect to plane 1. 
Nodes lying on plane 1 are fixed in the X direction. Plane 3 is fixed in Z direction and 
the load is applied on plane 2. Rigid bcdy motion in Y direction is restricted by fixing 
few nodes in the middle part of plane 3. 


4.3 RESULTS AND DISCUSSIONS 

As stated before, the model is analyzed for two different types of loading: tension and 
bending. Parametric study is' carried out for the two cases. The effect of thickness of weld 
is also studied. 

For the parametric study, the semi major axis a of the elliptical crack is fixed at 
20 mm and the value of aspect ratio b/a is varied from 0.3 to 0.8. The /-integral is 
evaluated over two different domains with 2 and 3 layers of elements around crack front 
for each case. 


4.3.1 Results for tension loading 

The results are obtained for tension loading for the semi elliptical cracks of aspect ratio 
b/a equal to 0.3, 0.4, 0.6 and 0.8 are plotted. The analytical relations for K, for infinite 
plate having semi-elliptical edge crack is given by 


K,= 


1.12o- 0r&) 1/2 


sin 2 ^+ — cos 2 ^ 

\a) 


1/4 


(4.1) 


where, 




\ 1/2 


sin 2 <j> 


(4.2) 


and <j) is the angle made by the line joining the node on ellipse to the centre of ellipse with 


X axis[37]. 

Since the diameter of the pipe (267 mm) is large compared to the crack length 
(semi major axis a =20 mm ), a comparison can be made between the pipe results and the 
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analytical results for infinite plate. These comparisons are made in the following graphs. 
Graphs are plotted for the variation of ./-integral with the angle, where the angle 
represents the location of crack front node, and it is measured from the centre of the 
ellipse with respect to X axis in counterclockwise direction. 



Angle in degrees 


Figure 4.8 Variation of ./-integral along the crack front for b/ a = 0.3 and r = 133.5mm in tension 
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J-integral (in N/m) 23 J-integral (in N/m) 


1000 
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Figure 4.13 Variation of ./-integral along the crack front for b/ a — 0.8 and r = 25mm in tension 

It is clear from the graphs shown in Figs 4.8-4.13 that as the aspect ratio b/ a 
(and hence the crack depth, 6) increases, value of J-integral increases and also the 
variation of the J-integral along the crack front becomes nearly the same as that for the 
infinite plate. For b/a ratio less than 0.6, the value of J is maximum at the centre of the 

crack front (at 90°) and decreases as we move towards the free surface. But for b/a ratio 
equal to 0.6 or more, the J-integral value is comparable to the maximum value, even at 
the free surface of the pipe wall. So it can be concluded that for low b/a ratio (less than 
0.6), the crack has a tendency to grow in the radial direction which corresponds to the 
leak-before-break condition. As the b/a ratio increases (more than 0.6), the crack has a 
tendency of propagating circumferentially which corresponds to the full guillotine offset 
break (i.e., ins tantan eous circumferential fracture). So there must be a critical value of the 
aspect ratio up to which the leak-before-break condition prevails. 

Figures 4.12-4.13 are for a pipe of diameter 50 mm for two extreme aspect ratios 

(b/a = 0.3 and 0.8). From these plots it is clear that radius of pipe have a very little 
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influence on the value of ./-integral. For small bja ratio, the effect of diameter is 
negligible, while for larger b/a ratio it is slightly more. So practically, diameter has no 
effect on the J values when the thickness of the pipe and the aspect ratio b/a remain 
constant. 

4.2.1 Results for bending loading 

A similar parametric study is carried out for bending load with the same set of aspect 
ratios b/a and with the same value of semi major axis (a= 20mm) . The bending load is 
applied in a manner that it opens the crack. A bending moment of 25 kNm is applied 
about X axis on plane 2 in clockwise direction. The results are plotted in the following 
Figures. 



Figure 4.16 Variation of /-integral along the crack front for bja — 0.3 and r - 133.5mm in bending. 
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J -integral (in N/m) to J-integral (in N/m) 









Figure 4.17 Variation of ./-integral along the crack front for bj a = 0.8 and r = 133.5mm in bending 

It is clear from Figs. 4.14-4.17 that the trend is same as in case of tension and it is 
expected because the state of stress created by bending moment near the crack is same as 
that in tension. As the aspect ratio b/a increases, the ./-integral increases and the curve 
gets flatted (same as in the case of tension). 
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CHAPTER 5 

CONCLUSIONS AND SCOPE FOR FUTURE WORK 

5.1 CONCLUSIONS 

In the present work a Matlab code has been developed to evaluate the 2-D 7-integral. The 
results obtained for SENB, CT and CP specimen in plane strain condition have shown 
good agreement with analytical solutions. 

The 2-D code has been extended to evaluate the point-wise 7-integral along the 
curved crack front for the 3-D crack problems. The values of 7-integral of the 3-D crack 
problems have also been successfully validated with the analytical results pertaining to 
plane strain problems. 

Variation of point-wise 7-integral along the crack front in a cracked pipe having a 
semi elliptical crack has been obtained and parametric study lias been carried out by 
varying the aspect ratio and keeping the pipe geometry and semi major axis of the 
elliptical crack fixed. These parametric analyses have been carried out for two different 
types of loadings namely bending and tension. The trend obtained for tension case has 
also been compared successfully with for an infinite plate with semi elliptical edge crack 
in tension. The results obtained from bending analysis are having the same trend as in 
case of tension. 

5.2 SCOPE FOR FUTURE WORK 

The desired future extensions of this work have been listed below: 

1 . Currently the program uses elastic stress and strain to evaluate the 7-integral. The 
program can be extended to evaluate 7-integral for elasto-plastic large 
deformation problems by incorporating appropriate constitutive theories. 
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2. The effect of thermal strains can be included by introducing thermal strain terms 
in the /-integral expression. 

3. The /-integral for dynamic problems can also be evaluated by adding appropriate 
dynamic terms in the /-integral expression. 
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